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Abstract In this paper, we analyze the convergence of the preconditioned 
GMRES method for the first order finite element discretizations of the Helmholtz 
equation in media with losses. We consider a Laplace preconditioner, an in- 
exact Laplace preconditioner and a two-level preconditioner. Our analysis is 
based on bounding the field of values of the preconditioned system matrix in 
the complex plane. The analysis takes the non-normal nature of the linear 
system naturally into account and allows us to easily consider certain type of 
inexact Laplace preconditioners via a perturbation argument. For the two- level 
preconditioner, our convergence analysis takes into account a media, which has 
not been considered in previous works. 
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1 Introduction 

Efficiently solving the Helmholtz equation is one of the major challenges in 
the field of numerical analysis. Current methods have difficulties both with 
the discretization and with the solution of the resulting linear system when 
the frequency grows. The first of these two difficulties is due to the very large 
number of elements required to obtain a meaningful numerical approximation 
to an highly-oscillating function. For example, the analysis in [TSlfTO] states 



The author was supported by the grant 133174 from the Academy of Finland 



A. Hannukainen 

Aalto University, Department of Mathematics and Systems Analysis, P.O. Box 11100, FI- 

00076 Aalto, Finland 

Tel.: -1-358 9 470 23696 

Fax; +358 9 470 23016 

E-mail: antti.hannukainenCihut.fi 



2 



Antti Hannukainen 



that for the first order finite element method, the mesh size h should satisfy 
the bound k^/i <C 1 before the asymptotic convergence rate is achieved for a 
system with first order absorbing boundary conditions. Here k is the wavenum- 
ber related to frequency of the resolved problem. In practice, satisfying such 
requirement for large wavenumbers n leads to solving very large linear systems. 

Developing efficient solvers for the linear systems arising form discretiza- 
tions of the Helmholtz equation has proven to be considerably more difficult 
than for elliptic problems. This is due to the indefinite nature of the linear 
system. In addition, if absorbing boundary conditions or certain type of losses 
are included, the discretization matrix is also complex valued and non-normal. 
This and the indefinite nature of the problem render many successful solution 
methods for positive definite problems less useful for the Helmholtz equation. 
For example, when applied to linear systems arising from the Helmholtz equa- 
tion, multigrid methods suffer from problems both in the smoothing and in 
the coarse grid correction steps, sec e.g. [S]. To have a convergent method, the 
coarse grid has to satisfy the same density constraints as the original discretiza- 
tion. In practice this means that although working multigrid preconditioners 
have been developed, they only provide a small benefit over a direct solver. 

The current trend for solving the linear system is to use a preconditioner 
together with a suitable Krylov subspace solver. The indefinite matrix problem 
can be solved with several Krylov subspace methods, e.g. GMRES, BiCGStab, 
etc. (see |141[^ ). From these methods, a covergence theory exists only for the 
GMRES method (see e.g [U). Because of this, we consider different precon- 
ditioners in connection with the GMRES method. 

The existing preconditioners can be divided into two groups, shifted-Laplace 
preconditioners (see e.g [TPll^lSlfTT] ) and two- level preconditioners (see e.g. 
mis]). The shifted-Laplace preconditioners are further development of Laplace 
preconditioners, which have been studied by several authors, e.g. pTllTj. These 
methods arc successful in cutting the growth in the condition number due to 
the Laplace operator part. However, a K-dcpendcncy in the required number 
of iterations still remains for the preconditioned system. Based on numerical 
examples [T0ll9] , the number of iterations for the shifted-Laplace type precon- 
ditioners for problems with Dirichlet boundary conditions behaves like 0{k^) 
and like 0{k) for problems with absorbing type boundary conditions. Regard- 
less of this asymptotic behavior, the introduction of the shift-term leads to a 
lower number of iterations compared to pure Laplace preconditioners. As the 
iterative methods for solving indefinite systems are computationally costly and 
memory intensive, even a small reduction in the number of required iterations 
is important. 

The two-level preconditioners are based on combining a Laplace precondi- 
tioner with a coarse grid correction. These methods can deliver K-independent 
number of iterations, but they suffer from identical problems as multigrid 
methods. Namely, a direct solver has to be employed to compute the coarse 
grid correction on a mesh satisfying the same constraints with the original 
discretization. The analysis of such methods has been performed for real val- 
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ued Helmholtz equation in ^^SM by using tools from the analysis of additive 
Schwarz methods for elliptic problems. 

In this paper, we will develop a field of values (FOV) based method to 
analyze the convergence of the preconditioned GMRES for the finite element 
discretizations of the Helmholtz equation with homogeneous Dirichlet bound- 
ary conditions in lossy media. We will consider a Laplace preconditioner, an 
inexact Laplace preconditioner and a two-level preconditioner. A similar anal- 
ysis has been done in |12| for Hermitian positive definite split preconditioners 
using algebraic tools. The main difference to this work is that we estimate the 
FOV by using methods similar to the ones applied in the analysis of additive 
Schwarz preconditioners for elliptic problems (see e.g. pS]*). 

The convergence of GMRES with a Laplace or a shifted-Laplace precon- 
ditioner has been previously analyzed in the literature, e.g. [13], by using 
algebraic tools. For certain kinds of losses, the preconditioned system matrix 
is diagonalizable in the inner product induced by a weighted mass matrix. In 
this approach, the eigenvalues are analyzed and the non-normality is taken 
into account by considering the conditioning of the weighted mass matrix. 

Our analysis takes the non-normal nature of the problem automatically 
into account and allows us to analyze the inexact Laplace preconditioner via a 
perturbation argument. We are also able to give convergence bounds for two- 
level preconditioner in lossy media. As such a preconditioner is not positive 
definite nor Hermitian, it is not covered by previous works. We will also give 
a more detailed analysis of the K-dependency of the coarse grid mesh size H 
when using two-level preconditioners. Our analysis indicates that in the worst 
case, the constraint k^H <?C 1 should be satisfied to guarantee convergence 
of the two-level method. The same mesh size constraint is also valid for the 
actual computational grid. The different mesh size requirement in comparison 
to [TMT^ is due to different boundary conditions. 

The organization of this paper is the following. First we introduce the model 
problem and quickly review the field of values based convergence theory of the 
GMRES method. We then introduce the preconditioners and give bounds for 
the FOV in each case. We conclude the paper with numerical examples on all 
of the proposed methods. 

2 Preliminaries 

We consider the problem 



where k e M and polyhedral domain J? C M^, = 2, 3. For the analysis of the 
exact and inexact Laplace preconditioners we assume that a E M.,<j > 0. More 
general losses, a G L°°{n), 



Au + (k^ 



ia)u ~ f in S7 
u = on dQ 



(1) 



a > and a > a„i > in C ^2, 



(2) 
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are considered in connection with the the two-level preconditioner. The more 
general case allows the presence of lossless areas, but introduces additional 
challenges in the analysis. 

The weak form of problem ([1]) is: Find u £ Hl{fi) such that 

a(u, = yv£Hl[Q), (3) 

in which (•, •) is the L^(]7) - inner product and 

a{u, v) = (Vu, Vv) — K^{u, v) + i((Tu, v). (4) 

Under our assumptions on cr, the weak problem ([3]) has a unique solution. This 
follows from the unique continuation principle (see e.g. [U) and the Fredholm 
alternative. 

In the analysis of the two-level prcconditioners, we will use a duality argu- 
ment. For this purpose, we need to consider the regularity and stability of the 
solution to problem ([3]). The coarse grid mesh size requirement k^H will 
arise from the K-dependency of the stability estimate. The weak solution will 
have the same regularity as the Poisson problem. 

Theorem 1 Let f G L?'{Q), k G R, cr G L°°{n) and u he the weak solution 
to (QJ). Then u G iJ3/2+<5(|2) mth some (5 > 0. 

Proof Clearly u is also the weak solution to the problem 

Au = (k^ — i(j)u + / in i7, 
u = on dfl. 

This is a Poisson problem with the right hand side in L^{Q). Hence, the regu- 
larity of u follows directly from the regularity theory for the Poisson equation, 
see e.g. [T51I5]. 

The parameter 5 in the above theorem is dependent on the shape of the 
polyhedral domain f2. For example, in convex domains 5 = 1/2. This depen- 
dency is analyzed carefully in [T^[5] 

In addition to the above regularity result, we need the stability estimate 

\\uh/2+s < Csll/llo- 

Our interest is especially in the K-dependency of the constant Cs- Studying this 
dependency for a general a is very difficult. Hence, we will give the estimate 
only for the case cr G K. 

Theorem 2 Let f G L'^{fl), k G K, cr G M and let u he the weak solution to 
(QJ). Then there exist a constant C > 0, independent on n and a, such that 

IHl3/2+^<^^(l + ^) ll/llo. 

for some S e (0, 1/2]. 
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Proof Let {vili^i be the eigenfunctions of the Laplace operator, i.e., 

-Aipi = Xiipi, and {ipi,ipj) = Sij. 
In this basis, the solution to the Helmholtz equation is 



^-^ \i — + icr 



i=l 



The solution satisfies the Poisson problem 



^-^ Xi — K'^ + icr 



The L^(/2)-norm of the right hand side is 



Elementary computations give. 



+ (72 ^ 

< , > 0. 



By the spectral theory of Laplace opertor, Xi > 0,i > 0. Hence, we get. 



E(/'^o'<^i + - 11/11 
i=i ^ ^ 

Now, the regularity theory for the Poisson problem yields the desired estimate. 

In the following, we will consider solving the linear system arising from the 
finite element approximation of the weak problem ([3]) with first order elements, 
i.e., the finite element space 

Vh = {ue Hi {Q) I v\K e Pi [K) yK eTh}. (5) 

The triangulation or tetrahedralization Th is assumed to be quasi-uniform (see 
[2]). In this space, the finite element approximation is: Find Uh G Vh such that 

a{uh, v) = (/, v) \fv eVh- (6) 
This problem leads to the linear system 

Ax ^ b. (7) 

Under our assumptions on a and k, the matrix A e -^[W be non- normal, 

i.e.. 
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AA* ^ A* A. (8) 

and indefinite. 

In the following, we work with functions from the finite element space Vt 
and the corresponding coefficient vectors. All functions from the finite element 
space can be expressed as m = ^(x„)i<pi, where fi are the finite element basis 
functions. The vector of coefhcients for the finite element function u will be 
denoted by x„. Using this notation, the system matrix A is related to the 
sesquilinear form ^ via Aj,i = a((pi,ipj), i.e., we have 

a{u, v) = x*^x„ Vw, V eVh- 

In the following, the norm equivalence between the Euclidian norm of vec- 
tor x„ and the i^(]7)-norm of a function u is used. 

Lemma 1 Let u G Vh and f2 CZ M.'^. Then there exists positive constants 
c,C > 0, independent of h, such that 

c/i''|x„|' < < C/i'^|x„|' (9) 

Proof See, e.g., [2] 

3 Convergence of GMRES 

The GMRES algorithm [23 approximately solves the linear system 

^x b (10) 

by iteratively constructing the minimizcr of the residual, \Ax — b|, from the 
Krylov subspace IC„i = {h, Ah, A^h, . . . , A^^^^h}. In exact arithmetic, the 
GMRES algorithm finds the exact solution in at most n iterations for a general 
invertible matrix A £ C"^". 

The convergence of the GMRES iteration can be improved by using right, 
left, or split preconditioners. The split preconditioner is used, when the precon- 
ditioner matrix can be decomposed into two parts, e.g. by using the Cholcsky 
decomposition. As obtaining suitable decomposition for our preconditioner 
matrices is very costly, using split preconditioner is not an option in our case. 

Left preconditioning in GMRES leads to minimizing the residual 

\BAx~Bh\ . 

In our case, the preconditioner matrix B is ill-posed, so the iterative solution 
might be incorrect even for very small values of the residual. Due to this fact, 
we will only consider right preconditioning. In this case the residual to be 
minimized is 



\ABi-h\ 
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The actual solution is x = i3x, hence the right preconditioned GMRES mini- 
mizes the actual residual, regardless of the conditioning of the preconditioner 
matrix B. On a sufficiently fine grid, this quantity can be easily related to the 
error between the iterative solution and the exact solution to the linear system 
measured in the 7f^(l?)-norm. 

The convergence of the GMRES method is related to the minimization 
problem (see e.g. [Tl] ) 



in which is the residual on step i and Pi a monic polynomial of order i. 
As solving the minimization problem (jlip is far more costly than solving the 
original linear system, it is not a practical measure of the GMRES convergence 
rate. More useful bounds have been derived from pT|) in several alternative 
ways, depending on the properties of matrix A (see, e.g, A good 

comparison of different GMRES convergence criterions is given in [7| . 

If the matrix A is normal, i.e., it satisfies equation the convergence is 
characterized only by the location of the eigenvalues of A. When the matrix 
A is non-normal but diagonalizable, the eigenvectors of the matrix are not 
orthogonal and they have an effect on the convergence in addition to the 
location of the eigenvalues. The convergence of GMRES for general non- normal 
matrix equations can be also related to the properties of the pseudospectrum 
PP] or the field of values [H]. The FOV is defined as the set 



Due to the connection between vectors of coefficients and finite element func- 
tions, the FOV is naturally related to the properties of the sesquilinear form 
a{u,v). Hence, we have chosen to use a FOV based convergence criterion in 
our analysis. 

The convergence of GMRES is related to the dimensions and the location 
of the set (fT2|) in the complex plane. Several different convergence estimates 
based on the FOV can be derived. A simple estimate is given in [M], let 
D = {zgC| |z — c|<s}bea disc containing the FOV, but not the origin. 
Then, we have the convergence estimate 



In this work, we are mainly interested in studying the dependence of the con- 
vergence of preconditioned GMRES on the mesh size h as well as parameters 
a and k. As the bound p3p remains unchanged under scaling of the coordinate 
system, it will be sufficient to study how the relative size of the FOV depends 
on these parameters. 

The properties of the FOV have been extensively studied in the literature 
(see, e.g., [T71[H1ITF| ) . In the numerical results section, wc will compute FOV 



min |p(A)ro|, 



(11) 




(12) 




(13) 
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for the preconditioned linear systems by using the procedure from [T3]. This 
procedure is based on the the rotation property of the FOV, 

F{A) = e-'^^(e'^A) 

and on the fact that T{A) is located on the left half plane from the largest 
eigenvalue of 

H{A)^^{A + A*). 

By computing the largest eigenvalues for several rotated matrices e'^A, we will 
obtain a set containing J-'{A). 



4 Laplace preconditioner for a constant cr 

In this section, we consider using the finite element solution of the Poisson 
equation as a preconditioner for the Helmholtz equation. We will give bounds 
for the FOV of the preconditioned system, which gives a convergence estimate 
for GMRES via equation Our analysis is valid when the parameter a £ 
M,cr > 0. 

The preconditioner P : 14 — >■ V/j is defined as: For each u G 14 find Pu G 14 
such that 

(VFw, Vw) ^{u,v) V ve Vh- (14) 

The matrix form of the operator P is K^^M , where K is the stiffness matrix 
and M the mass matrix, i.e., 

x*i4rx„ = (Vu, Vi;) and x*A/x„ = (li, w) VM,uel4- (15) 
The right preconditioned linear system has the form 

AK-^M^ = b. (16) 

We immediately observe, that 

TilAK~^MyLu = a{Pu,u) MueVh- 

Using this connection, the FOV set of the preconditioned system can be written 
as 



F{AK-H4) = ['^^^^ x„eC", x„^oj 



(17) 



To give bounds for this set, we will study the sesquilinear form a{Pu, u) instead 
of working directly with the matrix AK~^M. A similar connection is also 
the basis for the derivation of the convergence estimates for additive Schwarz 
methods applied to elliptic problems. The main difference to our case is that 
we need to obtain estimates between the sesquilinear form and the Euclidian 



FOV analysis of preconditioners for the Helmholtz equation 



9 



vector norm. For elliptic problems, similar estimates are derived in the H^{Q)- 
norm. 

By the definition of the sesquilinear form (|4|) , we have 

a{Pu, u) = {\/Pu, Vm) - K^{Pu, u) + ia{Pu, u). (18) 
To bound the terms above, we use the following elementary result. 

Lemma 2 Let P he defined as in and u G Vh. Then 

{Pu,u)^\\VPu\\l (19) 
and there exist a constant C > 0, independent of h, a , and k, such that 

\\VPu\\o<C\\u\W 

Proof The equation (|19|) follows directly from the definition of the operator 
P, equation (fT4)) . 

\\VPu\\l = {u,Pu) ^ {Pu,u). 
Using the Cauehy- Schwartz inequality gives 

\\VPu\\l < \\u\\o\\Pu\\o. 
Applying the Poincare-Friedrichs inequality completes the proof. 

Using Lemma [21 the definition of P, and equation (fT8)) gives 

a{Pu, u) = \\u\\l - K^\\VPu\\l + ia\\VPu\\l. (20) 

Based on this equation, it is straightforward to derive bounds for the FOV set. 
We begin with the obvious bounds. 

Theorem 3 There exists a constant C > 0, independent of h, a, and k, such 
that 

T{AK-Hl) c [C(l - K'^)h'^, Ch''] X [0, Cah'']. 
where d is the spatial dimension. 

Proof By equation ((20)) . we have 

Rex^AK-^Myi^^ \\u\\l - k^\\V Pu\\l (21) 

and 

Im x^MAxu = cr|| VPm||^. (22) 
In addition, we clearly have 

ReK:AK-'x^^\\u\\l-^^VPu\\l<\\u\\l 
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and 

ImxlMAxu = cr||VPu||o > 0. 
Using Lcmma[2]to estimate ||VPu||o in (j2ip and (|22|) . gives 

RexlAK^^Mxu > (1 - Ck^)\\u\\1 

and 

/mx*MAx„ < Ccr||u||o. 
The result follows from combining the above estimates with Lemma [TJ 

These bounds state that the FOV set is located inside a rectangle that 
contains the origin. In such a case, the convergence estimate (jl3p does not 
deliver any information on the convergence of GMRES. Fortunately, we can 
improve the above bounds. 

Theorem 4 There exists a constant C > 0, independent of h, a and k, such 
that 

f ^2 2 ^ 

T{AK-Hl) C <^ z e C ch'^ Im z < Re z < Ch'^ Im z \ . (23) 

where d is the spatial dimension. 

Proof For each z E J'{AK~^M) there exists a x„ <E C" and a corresponding 
u GVh such that 

^ ^ -KlAK-^Mxu ^ a{Pu,u) ^^4) 

Taking the imaginary part of the above and using (|20p yields 

\\\/Pu\\l - :^x:x„. (25) 
Combining this equation with the real part of (|24p and using (|20p . we obtain 

^^^^Rea{Pu,u) ^^_^,Irnz ^^6) 

Applying the Lemma [T] completes the proof. 

Theorems [3] and [3] state that the FOV set for the preconditioned system is 
located at the intersection of a strip and a rectangle. The intersection does not 
contain the origin, hence it can be used in connection with to give conver- 
gence estimates for the preconditioned GMRES method. As all the dimensions 
of the FOV have an h'^ - dependence, the relative size of the FOV does not 
change when the mesh is refined. Thus, GMRES for the preconditioned sys- 
tem will converge with the same rate independently of h. The convergence rate 
depends on the distance of the FOV set from origin and the size of the set. 
These parameters depend on k and cr, which will determine the conver- 
gence speed. 
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5 Inexact Laplace Inverse 

In this section, we consider replacing the exact solution of the linear system 
arising from discretization of the Poisson problem in the preconditioner (jl6p 
by an approximate solution. The presented analysis is based on a perturbation 
argument and it is valid for inexactly solving (jl4p using a symmetric iterative 
method convergent in the || ■ ||o and || ■ |ji norms. Example of a method fitting 
to this category is the multigrid (MG) method (for simple convergence proofs, 
see e.g. [2]). Our analysis indicates that the Poisson problem should be solved 
more accurately for large values of the parameter k for guaranteed convergence 
of the GMRES method. 

In the following, we are solving the linear system 

ATx = b (27) 

by using an iterative method. We shall denote one iteration cycle as and 
N cycles as . All iterations start from a zero initial guess, so that the 
error after N steps is e^r = x — K~^h. Equation (P7|) gives 

gn = {k-^ - K-^)Kx. 

This motivates us to define an error propagation operator En :Vh^Vh such 
that 

Cat = Emu. 
The matrix form of this operator is 

{R-'^ - K-^)K. (28) 
In the following, we assume that there exists constants 70 and 71, 

as well as a constant C > 0, independent on 70 and 71, such that 

\\Enu\\i<C-i^\\u\\i and \\Enu\\q < C^i^ \\u\\^ -iu^Vh. (29) 

This simply means that the applied iteration converges in the H^{Q)- and 
I/^(J?)-norms. Such an assumption is directly satisfied by several iterative 
methods, e.g., by the multigrid method. 

We will denote the inexact preconditioner as P. The matrix form of this 
operator is 

The FOV set for the preconditioned system satisfies 
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We estimate the size of this set by combining the resuhs from the previous 
section with a bound for the perturbation set 

T{A{K-^ - K-')M). (30) 

Using operator notation, bounding the perturbation set translates to giving 
bounds for 



a{iP-P)u,u) = 

(V(P - P)u, Vu) - k{{P - P)u, u) + ia{{P - P)u, u). (31) 

The last two terms in the above equation are estimated with the following 
lemma. 

Lemma 3 There exists a positive constant C > 0, independent on k,<7, 70 
and 71 , such that 

\{{P~P)u,u)\<Cj^\\u\\l VueVh. 
Proof Using the Cauchy-Schwartz and the Poincare-Friedrichs inequalities gives 

\i{P~P)u,u)\<C\\V{P-P)u\\o\\u\\o. 
The semi-norm above can be evaluated as 

l|v(P-PM.= »pM£_2^. 

Using the matrix form of the error propagation operator (|28p and symmetry 
as well as definitions ([TB| and ([5]) gives 

(V(P - P)u, Vv) = ^IK{K-^ - k-^)Mxu = (it, Env). (32) 

Hence, 

||V(P-P).||„=sup(^. 
By the Poincare-Friedrichs inequality and the assumption (|29p 

||V(P - P)u\\o = sup < C7f hllo, 

vev^ ||Vu||o 

which completes the proof. 

Estimating the first term in the right hand side of pip is straightforward. 
The matrix form of the error propagation operator (j28p and Cauchy-Schwartz 
inequality yield 

(V(P - P)u, Vu) = («, Eu) < \\Eu\\o\\u\\o < l^\\u\\l (33) 
Combining this equation with Lemma |3] yields the following theorem. 
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Theorem 5 There exists a constant C > 0, independent 0/70, 71, k, h, and 
a, such that 

Re J^{A{K-^ - K-')M) C [-Ch^i-f^ + ^^^f ), C;i'^(7o^ + )] 

and 

Im HAik-"" - K-^)M) C [-Ch\^^ + a^^),Ch\^^ + ^7^ )] . 

where d is the spatial dimension and N the number of iterations. 

Proof We will use the connection between the FOV set and the sesquilinear 
form to derive the bound. First, we observe that 

{{P - P)u,u) = xlAI{k-^ - K-^)Mxu 
as the applied iterative method is assumed symmetric, the term 

is real valued. Hence, 

Re T{A{k-^ - K-^)M) = Re (V(P - P)u, Vu) - k^{{P - P)u, u) 

and 

Im F{A{K-^ - K-^)M) = Im (V(P - P)u, Vu) + cr((P - P)u, u). 

Combining these equations with the equation ([33]), Lemma [3] and Lemma [1] 
completes the proof. 

Theorem [5] states that the perturbation set is located inside a rectangle. 
The dimensions of this rectangle are dependent on the number of iterations 
taken with the iterative scheme and on the parameters a and n. The estimate 
states, how the size of the perturbation set F{A{K~^ — K~^)M) converges 
to zero when the number of iterations N is increased. 

From theoretical point of view, the implication of Theorem [S] is that the 
number of iterations should be increased when the parameter k grows to keep 
the size of the perturbation set small and the origin outside the actual FOV. 
In section 7, we will numerically demonstrate that the presented analysis cap- 
tures the behavior of the perturbation set. The number of GMRES iterations 
required to solve the preconditioned problem is strongly affected by the num- 
ber of iterations N . However, the inclusion of the origin into the FOV seems to 
have a very small effect on the convergence of GMRES. This reflects the fact 
that the convergence estimate (|13p is only an upper bound for the required 
number of iterations. 
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6 Two-level preconditioner for a general <t 

In this section, we consider a two-level preconditioner combining a solution 
to the Poisson problem with a coarse grid correction. Our analysis allows the 
parameter cr to be a function satisfying the assumption given in equation ([2]) . 

The idea behind the two-level preconditioners for the Hclmholtz equation 
is to use a coarse grid solution to eliminate all eigenfrequencies corresponding 
to an eigenvalue with a negative real part. Then a Poisson problem is solved 
to approximate the solution components corresponding to eigenvalues with a 
positive real part. Such methods are successful in keeping the number of itera- 
tions constant, but impose same mesh size constraints as for the computationl 
grid also on the coarse grid. 

In previous works on this type of preconditioners P^lll] , the K-dependcncy 
of the coarse grid mesh size has not been explicitly studied. A widely acknowl- 
edged rule of thumb is that the coarse grid has to satisfy the same mesh size 
requirements as the original grid. We will verify this rule by giving bounds for 
the coarse grid mesh size H depending on k, a and the stability constant Cs- 
By Theorem[51 the stability constant for cr G K is Cs ~ n^cr^^ . As we will see, 
in the worst case this estimate leads to the requirement 

tC'H < 1 

for the coarse grid mesh size H . Based on our analysis, it is easy to see that the 
same mesh size constraint should also be valid for the actual computational 
grid. One should note that the bound k^/i <C 1 given in P^ITO] is obtained 
for a problem with different boundary conditions and hence a different sta- 
bility estimate. In addition, our analytical results consider only the so-called 
asymptotic range, whereas the more elaborate analysis in [TSlfTO] take also into 
account the pre- asymptotic range, i.e. the case when the mesh size requirement 
is not satisfied. 

Our preconditioner has the matrix form 

RhA^^rJjM + K-\I - ARHA-^Rfi)M, (34) 

where Ah is the system matrix from the space Vh and K as well as M are 
as defined in ([TC)). The matrix Rh is the prolongation operator from Vh to 
Vh- The mass matrix is included to the preconditioner (p4)) to aid in for- 
mulating it as an operator in the following analysis. The term RhAJj^RJjM 
corresponds to solving the original problem in the space Vh and the term 
K~^{I — ARhAJj^RJj)AI to solving the Poisson problem for the residual. 

As in the previous sections, we will derive an estimate for the FOV of the 
preconditioned linear system. In order to do this we will first interpret the 
preconditioner ((M)) as a sum of two operators, Ph and Q. The operator Ph 
corresponds to the solution of the original problem in the coarse space and Q 
to the solution of the Poisson problem for the residual. 

The operator Ph ■ Vh ^ Vh is defined as: For each u & Vh find Phu G Vh 
such that 
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a{PHU,VH) = {u,vh) ^vh^Vh- (35) 

The matrix form of this operator is RhA^R^M. To define the operator 
Q, we first introduce an a{u, v) - orthogonal projection operator Ih ■ Vh — > Vh- 
This operator is defined as: For each u £ Vh find Ihu € Vh such that 

a{v,lHu) ^ a{v,u) \/veVh- (36) 

The matrix form of this operator is RhAJj*RJjA*. An immediate consequence 
of this definition is the orthogonahty property 



a{vH,{I-lH)u) = VveVn. (37) 

This property will be used frequently in the following analysis. 

The term K~^{I — ARhA'^ R^)M in the prcconditioner (p4|) corresponds 
to the operator Q Vh ^ Vh defined as: For each u £ Vh find Q £ Vh such 
that 



{VQu,^v)^{u,{I - Ih)v) Vt;el4. (38) 

The main task in analyzing the two-level preconditioned is to derive a con- 
vergence result for the operator Q. This result is obtained by studying the 
properties of the projection operator Ih- The applied techniques are rather 
standard for the analysis of the finite element discretizations of the Helmholtz 
equation. However, as we want to explicitly state the dependency of the re- 
quired coarse grid mesh size on if, k, cr, and Cs some modifications to the 
approach used in previous works, ^SlI^ H] has been made. In the derivation 
of the convergence estimate, we will use the standard nodal interpolation op- 
erator to the coarse space, -kh- For a function u e H^/"^^^ {Q). 5 S (0, 1/2] this 
operator has the approximation property: 

\\u-nHu\\i<CH^/^+^\u\:i,2+s and Wu-nnuh < CH''/^+^\u\^,2+s, m 

where the constant C > is independent of H . For a function u e H^{S1), the 
convergence is obtained only in the L^(^?)-norm 

\\u-TrHu\\o<CH\\Vu\\o, (40) 

where the constant C > is independent of the mesh size H. Proofs of these 
approximation results can be found, e.g., from [2]. We begin with a standard 
convergence estimate for Ih- 

Lemma 4 There exists a constant C > 0, independent of h, H , Cs, k and 
cr,0 such that 

for some S G (0, 1/2] and H sufficiently small. 



16 



Antti Hannukainen 



Proof The proof is based on the standard duahty argument (see e.g. [5]). The 
L^-norm can be evaluated as 

\\(I - Ih)u\\o = sup — . (42) 

vevn \m\o 

The dual problem is: Find (p G Hq{Q) such that 

a{ip,w) = {w,v) yweH^if2). (43) 

Choosing the test function as w = (I — Ih)u and using the orthogonality 
property p7p gives 

{v, (/ - Ih)u) = a{ip - TTHip, (/ " Ih)u), 

where tth^^ G Vh is the nodal interpolant of ip. The interpolation error esti- 
mates given in equation (|39p and regularity Theorem [T] gives 



Hip ~ ^p, (I ~ In)u)\ < CCsH^/^+'\\V{I - Ih)u\\o\\v\\o 

+ CCsin' + ||a||^)ij3/2+^ 11 (/ _ Ih)u\\o\\v\\o. 

Combining this estimate with (I42p and reorganizing the terms yields 

II (/ - IhMo < ,^cCs{n^ + \\a\UHV^+^^^^^' 

which is valid for H such that 1 - CCs{k^ + ||cr||oo)i^^^/^+'' > 0. 

The above Theorem gives a convergence result for the operator Ih if the 
coarse grid satisfies the condition 

The parameter S in the estimate depends on the regularity of the solution to 
the dual problem As the load for the dual problem is always from the 

space H^{f2), the regularity depends completely on the shape of the domain. 
By Theorem [1] this dependency is the same as for the Poisson problem. In 
the following, we will state the coarse gird mesh size requirements under the 
assumption ||tT||oo ^ Under this assumption, the iJ^(J7) regularity of the 
dual problem, and in the light of the stability estimate given in Theorem ^ 
the mesh size requirement given in Lemma |4] translates to 

K^H < 1. 

In the following, we will use Lemma 2] in a shorter form: There exists a 
positive constant C > such that 

II (/ - IhHo < CCsH^'^^'\\V(I - Ih)u\W 
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The first implication of Lemma E] is tliat tlic scsquilinear form 

aiiI^lH)u,iI^lH)u) 

behaves in some sence as an coersive operator. This resuhs is required in the 
following analysis to obtain a K-independent bound for the term || V(/— ///)u||o- 

Lemma 5 Let u G H^{Q) and let the coarse grid mesh size H he such that 

K^CsH^^^+^ < 1 and n^ClH^+^^ < 1. (44) 
Then there exist a constant a > 0, independent on h,H ,Cs ,k, and a, such that 

Re - Ih)u, (/ - Ih)u) > a|| V(/ - Ih)u\\1 
Proof By the definition of the sesquilinear form, we have 

Re a{{I - Ih)u, (/ - Ih)u) = || V(/ - Ih)u\\1 - k^UI - Ih)u\\1 
Due to the assumptions (HH), we can apply Lemma H) This yields 

Re a((/ - Ih)u, (/ - Ih)u) > (1 - CK'ClH'+^')\\ViI - Ih)\\1 
Hence, the coarse grid mesh size has to satisfy the requirement 

1 - Ck^CIH^+^^ > 0. 

This is, K^cliji+a* « 1. 

The major implication of the above result is the requirement 

For a € R,a <^ and the H'^[f2) regularity of the dual problem (^5)) . this 
implies the requirement K^a~^H^ ^ 1. Hence, in the worst case the coarse 
grid has to satisfy the requirement 

K^iJ < 1. 

A second important consequence of Lemma [5] is the k- independent bounded- 
ness of the operator lu- This result is required, when deriving convergence 
estimates for the operator Q. 

Lemma 6 Let the coarse grid mesh size H be such that the assumptions |^^[ ) 
are satisfied and let u G H}^{f2). Then there exist a constant C > 0, indepen- 
dent on h,H ,Cs, and a, such that 



llv(j-/ff)w!io <q|Vu||o. 
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Proof Under the assumptions (jH]) , Lcmnia[5]states that there exists a constant 
a > 0, independent of h, H, Cs, and a, such that 

all V(/ - Ih)u\\1 < Re a{{I - Ih)u, {I - Ih)u). 

Now, let TTff It G Vh be the nodal interpolant of u. By the orthogonality prop- 
erty given in equation p7|) . we have 

q|| V(/ - Ih)u\\1 < Re a{{I - t:h)u, (/ - Ih)u). 

That is 

a\\V{I-lH)u\\l < Re {V {I - 7Th)u,V {I - Ih)u) 

- tih)u, {I - Ih)u) + i{<j{I - tih)u, {I - Ih)u)) . 

Using the Cauchy-Schwartz inequality, the interpolation error estimate (j40p . 
and the boundedness of the interpolation operator, we get 

a\\V{I-lH)u\\l < ||Vu||o||V(/-/ff)u||o+C(||a||,oW)i/||Vu||o||(/-/ff)ii||o. 

Applying the Poincare-Friedrichs inequality and dividing by ||V(/ — Ih)u\\o 
completes the proof. 

A direct consequence of Lemmas |4] and [6] is an approximation property for 
the operator Q. The approximation property is obtained both in the H-^{f2)- 
and i^(/2)-norms. The convergence result in the i^(]7)-norm follows from the 
H^{f2) approximation property and the duality argument. 

Lemma 7 Let the coarse grid mesh size H be such that the assumptions given 
in equation are satisfied. In addition let u € Vh- Then there exist a con- 
stant C > 0, independent of h,H ,Cs , k, and a, such that 

llVg^llo < CCsi?^/'+*||w||o yueVh. (45) 
for some S G (0, 1/2]. 

Proof By the definition of operator Q, equation we have 

\\VQu\\1 = {u,{I-Ih)Qu). 
The Cauchy-Schwartz inequality and Lemma S] yield 

\{u, {I - Ih)Qu)\ < CCsH^^''+'\\u\\o\\W{I - lH)Quh, 
Applying the boundedness result from Lemma |6] completes the proof. 
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Lemma 8 Let the coarse grid mesh size H be such that assumptions given in 
equation ji44\ l are satisfied. In addition let u € Vh. Then there exist a constant 
C > 0, independent of h,H ,Cs, n, and a, such that 

IIQwIlo < CC5-ffi+'*||u||o "^u^Vh. (46) 
for some 6 £ (0, 1/2]. 

Proof Let ip G HQ{f2) be such that 

As Qu € Vh, also Qu € L^{fl). Hence, by eUiptic regularity theory, € 
iy3/2+'5(/2) and ||<(a||3/2+5 < C||Qu|jo. By the definition of Q, we have the 
orthogonality property 

(V7rff(^, VQu) = 0. 
By this property and definition of opeartor Q, equation (|38p . we have 

(V((/3 - tthv?), VQu) ^ IIQujIo- 
Next, the Cauchy-Schwartz inequality leads to 

||Q^^||^< llV(^-7r^^^)||o||VQ«||o 
Applying the interpolation error estimate given in equation p9p . the stability 
estimate for tp, and Lemma [7] gives 

\\Qu\\l< CCsH^+^'\\Qu\\o\\u\\o. 
Dividing with ||(5uj|o completes the proof. 

Now, we are finally at the position to give bounds for the FOV. The pre- 
conditioner ([M)) has the operator form Ph + Q^ so we need to bound 

a{{PH + Q)u,u). 

From the definition of the interpolation operator Ih given in equation pop and 
the definition of operator Ph given in equation psp . it immediately follows 
that 

a{PHU,u) = a{PHU, Ihu) ~ {u,Ihu). 
The definition of the operator Q given in equation yields 

a{Qu, u) ~ (u, (/ - Ih)u)q — K^iQu, u)o + [{oQu, u). 
Combining the two above identities gives 

a{{PH + Q)u, u) = ||u||o - K^iQu, u)o + \{aQu, u). (47) 

Estimating the last two terms above using Lemma |5] leads to the desired 
bounds for the FOV. 
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Theorem 6 Let the coarse grid mesh size Hsatisfy the assumptions 

k^ClH^+^^ < 1 and k^CsH^^^+^ < 1. (48) 
Then there exists constants c, C > 0, independent of h, H , k, and a , such that 

ReFc [ch" (1 - Cs{k'' + \\a\\^)H^+^') , Ch" (l + CsiK^ + \\a\U)H^+^')] 
and 



where d is the spatial dimension and 5 G (0, 1/2]. 
Proof By equation (|T7)). we have 

Re a{{PH + Q)u, u) = \\u\\l + Re (i(crQu, u) - k^{Qu, u)) (49) 

and 

Im a{{PH + Q)u, u) = Im (i{aQu, u) — k^{Qu, u)) . (50) 

Under the assumptions made in equation (|48)). we can use the convergence 
estimates for term ||Qu||o given in Lemma |8] to bound the above terms. This 
leads to the estimates 

\{Qu,u)\<CCsH^+^'\\u\\l 

and 

\i<jQu,u)\ < \\<j\\^CCsH^+^'\\u\\l 

Combining these inequalities and Lemma [1] with equations (|49|) . ((50)) . com- 
pletes the proof. 

The bounds given in the above theorem reflect the requirement on the 
approximation properties of the space Vh ■ The coarse grid has to be sufficiently 
dense, before the FOV set is completely located at the right-half plane. Based 
on the above result, this happens in a convex domain when k^H <^ 1. However, 
to obtain the estimates applied to derive the result the coarse grid mesh size 
requirement k^H <^ 1 was made. Hence, the coarse grid mesh should satisfy 
the very strict requirement k^H ^ 1. 
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Fig. 1 The first two mesh levels for the two dimensional test ease. 



7 Numerical Examples 

In this section, we study the presented theory by numerical examples. We 
verify the bounds derived for the FOV by computing the actual sets with the 
procedure presented in Section 3. We will also solve the preconditioned linear 
system to show the relationship between the FOV and the number of GMRES 
iterations required to solve the problem. 

The numerical tests arc composed of two parts. First, we will consider the 
two dimensional unit square, f2 = (0,1)^. In this test, we use a family of 
triangulations composed by uniformly refining an initial mesh. The first two 
mesh levels are presented in Fig. [T] 

The second part is a computationally more realistic three dimensional test 
case. As the number of degrees of freedom required in a three dimensional 
domain is considerably larger compared to the two dimensional case, we have 
not computed the FOV. The main focus is in the convergence of the precon- 
ditionerd GMRES method. 



7.1 Laplace preconditioner 

We begin by studying the exact Laplace preconditioner presented in Section 
4. Our aim is to verify the bounds given for the FOV in Theorems [3] and E) 
These theorems state that the FOV for the Laplace preconditioned system in 
two dimensions is contained inside a rectangle (— ck^/i^, C/i^) x (0, Cah^) and 
the strip 

( ^2 ^2 > 

zeC ch^ Imz<Rez<Ch^ Im z \ , (51) 

(7 ) 

in which constants c, C > are independent of k, a and h. 

The ft,^-scaled FOV for the Laplace preconditioned system are presented 
in Fig. [2] for (7 = 1 and three diff'erent parameter values = 1, 10, 50. Each 
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FOV set for Laplace preconditioner, ^1 ,o=1 FOV set for Laplaoe preconditioner, =10,o=1 




0.1 0.3 0.5 0.1 0.3 0.5 

FOV set for Laplace preconditioner, =50,o=1 



0.03 



0.02 



0.01 








-1 -0.5 0.5 

Fig. 2 The h^-scalcd FOV for the exact Laplace preconditioner in the two dimensional test 
case. The sets arc computed on mesh levels from one to five. The parameters are cr = 1 and 
= 1,10,50. 

figure contains the FOV from mesfi levels one to five. Based on these results, 
the /i^-scaled shape of the FOV seems to be dependent only on the parameters 
K and (T. 

To verify the theoretical bounds for the rectangle bounding the FOV set, 
we have computed the /i^ - scaled dimensions of this rectangle for mesh levels 
from one to four and for different parameter values. The bounds are presented 
in Figs. 13] and One can immediately observe from these results that the 
dimensions of the rectangle converge to a limit value when the mesh level 
is increased. The conclusion is that the /i^-scaled size of the rectangle stays 
constant, as predicted by the analysis. Based on Figs.[3]and|4l the lower bound 
of the real part depends linearly on and the upper bound of the imaginary 
part linearly on a. The upper bound for the real part as well as the bounds 
for the imaginary part are -independent. The bounds for the real part are 
cr-independent and the lower bound for the imaginary part is very close to 
zero. These results are in good accordance with Theorem [31 

Next, we consider the strip containing the FOV for different parameter 
values and mesh levels from one to four. The strip is computed by finding the 
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Laplace prec, bounds for Im part of FOV set, o=5 
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Fig. 3 The /i^-scaled bounds for the real and 
the FOV of the Laplace preeonditioned system. 



imaginary parts of the rectangle containing 
The parameter tr = 5. 
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Fig. 4 The /I'^-scaled bounds for the real and 
the FOV of the Laplace preconditioned system. 



imaginary parts of the rectangle containing 
The parameter = 1000. 



two lines with the slope ^ bounding the FOV from above and below. The 
bounds for the x-intercept of the strip ([51]) are visualized in Fig. [5] Based on 
these bounds, the h^-scaled x-intercept points of the two lines bounding the 
FOV are independent of k, cr and h. This is as predicted by Theorem U] 

To demonstrate the dependency of the required number of GMRES itera- 
tions for the Laplace preconditioned system on a and k, we have solved the 
problem with different parameter values and the load function / = 1. The 
level 7 mesh was used in the computations. The stopping criterion for the 
GMRES iteration was set to 10~^. The required number of GMRES iterations 
is visualized in Fig. [6] Based on these results, a linear dependency between 
the required number of iterations and is observed. The slope of the to 
number of iterations line is dependent on a. This is due to the distance of the 
FOV from the origin being dependent of the ratio of and a. 



24 



Antti Hannukainen 



Laplace prec, bounds for the strip, = 1 00 Laplace prec, bounds tor the sthp, - 500 



0.5 


S S: S: S: * 






"3 


0.5 


S S: 3: 3= 






= 3 


0.4 


•- -•- -• • 


- -•- 


• 


•- 


• 


0.4 






• 




• 


S "'^ 












S 0.3 












h -scaled s 


•- -• • • 

» 4= 4= 4 


•- 


• 

: = ^ = 


•- 

= =8= 


• 

= = 4 


'w 
■8 

S 0.2 
0.1 


» 4= 4= 4= 4= 






• 

= 4 




- • - Level 1 

- e - Level 2 

- ^ - Level 3 

- B - Level 4 












- • - Level 1 

- e - Level 2 

- ♦ - Level 3 

- B — Level 4 












1 5 10 15 


20 


30 

o 


40 


50 

Laplace prec. 


bounds for the strip, - 


1 5 10 15 
1000 


20 


30 

a 


40 


50 



»4=4=4=4= = = ^ = = ^ = = 4 



- ♦ - Level 1 

- e - Level 2 

- ♦ - Level 3 

- B — Level 4 

1 5 10 15 20 30 40 50 



Fig. 5 The h^-scalod bounds for x-intercepts of the strip containing the FOV for ex- 
act Laplace preconditioner in the first test case. The parameter k has the values = 
100, 500, 1000. 



0.4 
01 0.3 

■D 

S 0.2 



Laplace prec, number of GMRES iterations 



180 
160 




/ 


140 




/ « 


120 
100 




✓ ,11 

' i. 

✓ *, ^ . ■ 

■"I 


80 






60 


- • - <j= 10 - 


40 
20 




- e - o.ioo 

- ♦ - 0=200 - 

- B - o =300 
" J» ' ' 0=400 - 







■ ■< ■ ■ 0=500 



100 500 1000 1500 2000 



Fig. 6 The number of GMRES iterations required to solve the Laplace preconditioned 
system for / = 1. A seventh level mesh and stopping criterion 10"*^ were used in the test. 
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7.2 Inexact Laplace preconditioner 

Next, we will consider replacing the exact Laplace preconditioner with a multi- 
grid solver. The multigrid solver uses T^-cycle iterations with one pre- and 
postsmoothing step with a Richardson smoother. A mesh hierarchy of mesh 
levels from two to seven is used in the tests, if not otherwise stated. Our aim is 
to demonstrate the bounds for the FOV of the MG preconditioned system pre- 
sented in Section 5. The main interest lies in the behavior of the perturbation 
set 

when K, CT, the number of multigrid iterations, or the number of levels are 
varied. In Section 5, we have shown that the perturbation set is located inside 
a /i^-scaled rectangle 

(-c7o^ - c«^f , Cj^ + C.^f ) X (-c7o^ ^ ca^f , Cj^ + Ca'j^) 

where c, C > are constants independent of the computational mesh, k, a, 
and the applied iterative scheme. The parameters 70 and 71 are the error 
reduction factors for the multigrid solver and TV is the number of ^-cycles. 

The FOV for the MG preconditioned system is presented for parameters 

= 1, 10, 50, <J ~ 1, and a varying number of T^-cycles in Fig. [T] Based on 
this figure, one can verify that the FOV converges to a limit set, when the 
number of V^-cycles in increased. This behavior is as predicted in Section 5. 
Examples of the corresponding perturbation sets are given in Figure [SI The 
diameter of the perturbation sets clearly converges to zero when the number of 
F-cycles is increased. Due to the relatively small values of k, all perturbation 
sets look rather similar. 

Next, we will study the dependency of the dimensions of the perturbation 
set on the number of l/-cycles and h. In this test, the coarse mesh is kept 
unchanged, but the number of levels in the multigrid hierarchy is varied from 
two to seven. The coarse grid for each test is the level two mesh. The results 
are presented in Figs. 1^ and ITUl 

From Fig.[3J one can observe a linear dependency between and the upper 
bound for the real part of the perturbation set. No K^-dependency is observed 
for the imaginary part of the set. The number of levels in the multigrid hierar- 
chy affect both the slope and the y-intercept point of the computed lines. The 
bounds in Section 5 predict that the parameter 70 determines the y-intercept 
point and parameter 71 changes the slope of the lines. The error reduction 
factors 7o and 71 of the MG method depend on the number of levels, see e.g. 
[2], which explaines this phenomenon. 

The results for varying a in Fig. [TO] are very similar. The lower bound for 
the imaginary part of the perturbation set depends linearly on cr, whereas the 
real part is cr- independent. Again, the intercept points and slopes of the lines 
change due to the changing error reduction factors. 
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MG preconditioner, FOV set forK^=1 ,o=1 MG preconditioner, FOV set for K^=10,a=1 
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Fig. 7 The /i^-scaled FOV for the multigrid preconditioned system. Up to 8 multigrid 
V-cycles have been used. The parameter (7 = 1 and = 1, 10, 50. 



The results on the size of the perturbation set are in accordance with the 
derived bounds. However, the estimate for the upper bound of the real part 
and lower bound for the imaginary part seem to stay constant for varying 
K and a. The bounds given in Section 5 overestimate the size of the set by 
predicting k^- and a-dependency in these cases. As the FOV is a direct sum 
of the perturbation set and the FOV set for the Laplace preconditioner, the 
derived bounds manage capture the behavior of the FOV regardless of this 
overestimation. 

Finally, we study how the dimensions of the perturbation set behave when 
the number of l/-cycles grows. We have used a hierarchy with four levels and 
varied the number of l^-cyclcs from one to four. The dependency of the bounds 
for the real part on n and the number of V^-cycles is visualized in Fig. [TT] and 
the dependency of the imaginary part on a and the number of T^-cycles in 
Fig. [121 In addition, we have visualized the development of the bounds for 
parameters = 1000 and cr = 5 for larger number of V^-cycles in Fig. [T31 

Based on these results, the size of the perturbation set decreases as pre- 
dicted when the number of l^-cycles is increased. The slope and the intercept 




points of the lines tend to zero. Based on Fig. [T31 the convergence speed is 
same for all dimensions of the set, as predicted by in Section 5. 

To study the dependency of the required number of MG preconditioned 
GMRES iterations on the parameters and the number of V^-cycles, we consider 
the problem with / = 1. The MG method uses a mesh hierarchy with six levels, 
with the level two mesh acting as the coarse grid. The stopping criterion for the 
GMRES iteration was set to 10"^. The required number of GMRES iterations 
for different parameter values and different number of \^-cycles is visualized 
in Fig. [H 

The interesting factor for the multigrid preconditioner is the dependency 
of the required number of GMRES iterations on the number of ^-cycles. This 
is demonstrated in Fig. [14] by comparing the number of iterations for varying 
number of multigrid V^-cycles against one using ten T^-cycles. Based on these 
results, one can observe a clear dependency between the number ^-cycles and 
the number of required iterations. For large values of k, the number ofV -cycles 
has a quite big effect to number of required number ofG MRES iterations. This 
is due to the K^-dependency of the size of the perturbation set. The inclusion 
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Fig. 9 The h-^-scaled bounds for the real and imaginary parts of the perturbation set for the 
MG preconditioner with a single V-cycle. A Level one mesh is used for the coarse grid and 
the number of levels in the multigrid hierarchy is varied from two to seven. The parameter 
is cr = 5. 
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MG prec, bounds for Re part of perturbation set, ic =1000 
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Fig. 10 The h'^-scaled bounds for the real and imaginary parts of the perturbation set 
for the MG preconditioner with a single V-cycle. A level two mesh is used for the coarse 
grid and the number of levels in the multigrid hierarchy is varied from two to seven. The 
parameter is = 1000. 



of the origin to the FOV for the MG preconditioned system seems to be quite 
irrelevant for the convergence. 
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Fig. 11 The h^-scaled bounds for the real and imaginary parts of the perturbation set for 
the MG preconditioner. The coarse grid and the number levels in the multigrid hierarchy 
are fixed and the number of V-cycles is varied. The parameter is cr = 5. 
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Fig. 12 The h^-scalcd bounds for the real and imaginary parts of the perturbation set for 
the MG preconditioner. The coarse grid and the number levels in the multigrid hierarchy 
are fixed and the number of V-cycles is varied. The parameter is = 1000. 
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MG prec, bounds for perturbation set, k = 1000, o=5 
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Fig. 13 Convergence of the h-^-scalcd bounds for the real and imaginary parts of the per- 
turbation set for the MG preconditioner as a function of the number of \/-cycles. The 
parameters are ct = 5 and = 1000 
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Fig. 14 The left figure shows the number of GMRES iterations required to solve the MG 
preconditioned system for / = 1. The right figure shows the convergence of the number of 
GMRES iterations by comparing number of iterations for varying number of l^-cycles to 
MG preconditioner using ten V-cycles. 
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Fig. 15 The h^-scaled FOV sets for the two level preconditioned system. A seventh level 
fine mesh was used. 

7.3 Two-level preconditioner 

Finally, wc consider the two-level preconditioner presented in Section 6. Our 
aim is to verify the bounds given in Theorem [S) The most interesting param- 
eter for the two-level preconditioner is the mesh size of the coarse grid. This 
parameter effectively determines the amount of computational work required 
to solve the linear system. 

In Fig. [ini the FOV is computed for parameters k = Att and cr = 7 by 
keeping the fine grid level fixed to seven and varying the coarse grid level. 
Based on these results, if a sufficiently small coarse grid mesh size is used the 
FOV is located at the right half-plane and does not contain the origin. In 
addition, the imaginary part of the FOV converges to zero. This behavior is 
as predicted in Section 6. 

The interesting result of Section 6 is the implication that the coarse grid 
mesh size in a convex domain should be such that the term k^H is small. 
We have computed the smallest real part of the FOV set for several different 
values of k using a level nine computational grid and different coarse grid 
levels. The coarse grid level required before FOV is located at the right half- 
plane is presented in the Fig. As only few datapoints were computed it 
is difficult to determine if the requiremenet for the coarse grid mesh size is 
necessary or not. Unfortunatelly the FOV sets are computationally expensive 
to find, so we are not able to give a conclusive example. 

In Fig. [T7] the behavior of the rectangle bounding the FOV of the two-level 
preconditioned system for parameters a = 7 and k = 47r, Gtt, IOtt is visualized 
as a function of the coarse grid mesh size H. Based on these results, after a 
sufficiently dense coarse mesh is reached the bounds for the imaginary part 
converge to zero. The bounds for the real part on the other hand converge to 
limit values. These results are in accordance with the behavior predicted for 
FOV in Section 6. 
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Fig. 16 The coarse mesh level required for the real part of the FOV set to be located in 
the right half plane. A level nine fine grid was used. 
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Fig. 17 The /i^-scaled FOV sets for the two level preconditioner. The parameter a = 1 and 
K = 4iT, 67r, IOtt. 
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Fig. 18 The number of GMRES iterations required to solve the two-level preeonditioned 
system. A level nine fine grid wras used and the coarse grid level was varied. 



To study the two-level preconditioned GMRES method, we solve the linear 
system using a level nine fine grid and varying coarse grid level. The load is 
/ = 1 and the stopping criterion is set to 10^^. The results are presented 
in Fig. [TSl These results are as predicted in Section 6. After a sufficiently 
large coarse grid level, the number of iterations stagnates to a limit value. The 
stagnation point depends on the parameter n. An interesting observation is 
that the number of iterations converges regardless of the inclusion of origin in 
the FOV. 



7.4 A three dimensional example 

We conclude the numerical examples by considering a three dimensional cube, 
fi = (0, 1)"^. A hierarchy of uniformly refined tetrahedral meshes with five levels 
has been used in all of the tests. The finest level mesh had approximately 
750 • 10'^ degrees of freedom and 4.5 • 10^ tetrahedral elements. As we have 
demonstrated the behavior of the FOV in detail in the two dimensional test 
case, our interest will be solely on the number of iterations required to solve 
the preconditioned systems using GMRES. The load for all three dimensional 
test cases is chosen as / = 1 and the stopping criterion for GMRES is chosen 
as 10-6. 

We begin our experiment with the exact Laplace preconditioner. To solve 
the Laplace problem, a preconditioned conjugate gradient (PCG) method with 
a single multigrid T^-cycle as a preconditioner was used. The stopping criterion 
for the PCG method is set to 10"^^. The parameters a and k are both varied. 

The required number of GMRES iterations for different parameter values is 
presented in Fig. [121 Based on these results, the number of iterations behaves 
qualitatively as predicted by the theory. The number of iterations is dependent 
both on and a. When a grows, the presented theory predicts that also the 
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Fig. 19 The number of GMRES iterations required to solve the Laplace preconditioned 
system in the three dimensional test case. 
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Fig. 20 The loft figure shows the number of GMRES iterations required to solve the 
three dimensional test problem. The right figure shows the convergence of the number of 
iterations by comparing the number of iterations for varying number of V-cycIes to the MG 
prcconditioner using ten V-cycles. 



distance from the origin to FOV increases. This can be seen as the changing 
slope of the K^to number of iterations Hnes in the Fig. [TOl 

Next, we replace the exact Laplace prcconditioner with a multigrid based 
prcconditioner. The interesting question here is the dependency of the number 
of GMRES iterations on the number of l/-cycles. The required number of 
GMRES iterations for different parameter values is presented in Fig. [20l Based 
on these results one can clearly observe that the required number of iterations 
is strongly dependent on the number of multigrid F-cycles. Again, the inclusion 
of the origin into the FOV is not observed in the required number of GMRES 
iterations. 

Finally, we consider the two-level prcconditioner. The mesh levels from 
one to four arc used as the coarse grid. The same procedure as for the exact 
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Coarse grid levels 


K 


1 


2 


3 


4 


4n 


52 


37 


28 


25 




197 


131 


73 


37 


IOtt 


n.c. 


n.c 


n.c. 


130 



Table 1 Number of GMRES iterations required to solve the problem using the two-level 
prcconditioner. The computational mesh has approximately 750-10'^ degrees of freedom and 
4.5 ■ 10® tetrahedral elements. The iteration is deemed not to converge, if more than two 
hundred iterations are required. 



Laplace preconditioner is used to solve the Poisson problem on the fine grid. 
The number of GMRES iterations for different values of k and different coarse 
grid levels are presented in Table [T] From the number of iterations one can 
qualitatively verify the result given in Section 6. Namely, large values of k 
require more dense coarse meshes before the number of iterations converges to 
a limit value. In the present case, this happened only for the parameter value 
K = 47r. Probably, even a more dense computational mesh would be required 
to resolve the solution for k = IOtt. 
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